#!/usr/bin/python
import numpy as np
eV = 1.6e-19
MeV = 1e6 * eV
GeV = 1e9 * eV
KeV = 1e3 * eV
ekmin = 50 * KeV
ekmax = 10.0 * MeV
alpha = 1.0 / 137
#me = 9.1e-31
me = 0.511 * MeV
r0 = 2.8179e-15
myz = 79.0
logekrange = np.mgrid[np.log(ekmin):np.log(ekmax):500j]
ekrange = np.exp(logekrange)
logxrange = np.mgrid[np.log(1e-4):np.log(0.999):500j]
x_range = np.exp(logxrange)
np.savetxt("ele_ek.txt", ekrange)
np.savetxt("xrange.txt", x_range)
def fz(Z):
    z = alpha * Z * alpha * Z
    return 1.202 * z - 1.0369 * z * z + 1.008 * z * z * z / (1.0 + z)

def phi1(gamma):
    return 20.863 - 2.0 * np.log(1.0 + (0.55846 * gamma)**2) - 4.0 * (1.0 - 0.6 * \
            np.exp(-0.9 * gamma) - 0.4 * np.exp(-1.5 * gamma))
def phi2(gamma):
    return phi1(gamma) - 2.0 / 3.0 / (1.0 + 0.65 * gamma + 6.0 * gamma * gamma)
def psi1(epsilon):
    return 28.340 - 2.0 * np.log(1.0 + (3.621 * epsilon)**2) - 4.0 * \
            (1.0 - 0.7 * np.exp(-8.0 * epsilon) - 0.3 * np.exp(-29.2 * epsilon))
def psi2(epsilon):
    return psi1(epsilon) - 2.0 / 3.0 / (1.0 + 40.0 * epsilon + 400.0 * epsilon**2)

def lrad(Z):
    return np.log(184.1 / Z**(0.3333))

def lradp(Z):
    return np.log(1194.0 / Z**(0.6666))

def cs(Z, x, ek):
    delta = me * me / (2.0 * ek * x * (1.0 - x))
    gamma = 200.0 * delta / (me * Z**0.3333)
    epsilon = 200.0 * delta / (me * Z**0.6666)
    return alpha * r0 * r0 * ((4.0 / 3 * x * x - 4.0 / 3 * x + 1.0) \
            * (Z * Z * (phi1(gamma) - 4.0 / 3 * np.log(Z) - 4.0 * fz(Z)) + \
            Z * (psi1(epsilon) - 8.0 / 3 * np.log(Z))) - 2.0 / 3 * x * (1.0 - x)\
            * (Z * Z * (phi1(gamma) - phi2(gamma)) + Z * (psi1(epsilon) - psi2(epsilon))))
def cs_br(Z, y, ek):
    return 4.0 * alpha * r0 * r0 * ((4.0 / 3 - 4.0 * y / 3 + y * y) * \
            (Z * Z * (lrad(Z) - fz(Z)) + Z * lradp(Z)) + 1.0 / 9 * (1.0 - y) * \
            (Z * Z + Z)) / y
dx = x_range[1] - x_range[0]
sigma = np.zeros((ekrange.shape[0], x_range.shape[0]))
sigma_max = np.zeros(ekrange.shape)
def generatez(myz):
    for i in np.arange(ekrange.shape[0]):
        for j in np.arange(x_range.shape[0] - 1):
            dx = x_range[j+1] - x_range[j]
            sigma[i][j + 1] = cs_br(myz, x_range[j], ekrange[i]) * dx + sigma[i][j]
    np.savetxt("csbr" + str(int(myz)) + ".txt", sigma)

for zzz in np.arange(79, 80, 1):
    print "generating ", zzz
    generatez(zzz)

